/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of UHM (Unassembled HyperMatrix) sparse direct solver.
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#include "uhm/common.hxx"
#include "uhm/const.hxx"

#include "uhm/matrix/uhm/matrix.hxx"
#include "uhm/matrix/uhm/fla.hxx"

namespace uhm {
  extern double flop;

  static int solve_qr_1_flat( int fs, int ss,
                              linal::Flat_ ATL, linal::Flat_ ABL, 
                              linal::Flat_ t,   linal::Flat_ b,
                              linal::Flat_ T );


  static int solve_qr_2_flat( int fs, int ss,
                              linal::Flat_ ATL, linal::Flat_ ATR, 
                              linal::Flat_ t,   linal::Flat_ b );
  

  static int solve_qr_1_hier( int fs, int ss,
                              linal::Hier_ ATL, linal::Hier_ ABL, 
                              linal::Hier_ t,   linal::Hier_ b,
                              linal::Hier_ T );


  static int solve_qr_2_hier( int fs, int ss,
                              linal::Hier_ ATL, linal::Hier_ ATR, 
                              linal::Hier_ t,   linal::Hier_ b );
  
  
  void Matrix_FLA_::solve_qr_1_x() {

#ifdef UHM_HIER_MATRIX_ENABLE
    // ----------------------------------------------------------
    // ** Hier-Matrix 
    // ----------------------------------------------------------
    solve_qr_1_hier( this->fs, this->ss,
                     this->hier.ATL, this->hier.ABL,
                     this->hier.xt,  this->hier.xb,
                     this->hier.T );

#else
    // ----------------------------------------------------------
    // ** Flat-Matrix 
    // ----------------------------------------------------------
    solve_qr_1_flat( this->fs, this->ss,
                     this->flat.ATL, this->flat.ABL,
                     this->flat.xt,  this->flat.xb,
                     this->flat.T );


#endif
  }

  void Matrix_FLA_::solve_qr_2_x() {

#ifdef UHM_HIER_MATRIX_ENABLE
    // ----------------------------------------------------------
    // ** Hier-Matrix
    // ----------------------------------------------------------
    solve_qr_2_hier( this->fs, this->ss,
                     this->hier.ATL, this->hier.ATR,
                     this->hier.xt,  this->hier.xb );
#else
    // ----------------------------------------------------------
    // ** Flat-Matrix 
    // ----------------------------------------------------------
    solve_qr_2_flat( this->fs, this->ss,
                     this->flat.ATL, this->flat.ATR,
                     this->flat.xt,  this->flat.xb );

#endif
  }    

  void Matrix_FLA_::solve_qr_1_r() {
#ifdef UHM_HIER_MATRIX_ENABLE
    // ----------------------------------------------------------
    // ** Hier-Matrix 
    // ----------------------------------------------------------
    solve_qr_1_hier( this->fs, this->ss,
                     this->hier.ATL, this->hier.ABL,
                     this->hier.rt,  this->hier.rb, 
                     this->hier.T );
#else
    // ----------------------------------------------------------
    // ** Flat-Matrix 
    // ----------------------------------------------------------
    solve_qr_1_flat( this->fs, this->ss,
                     this->flat.ATL, this->flat.ABL,
                     this->flat.rt,  this->flat.rb, 
                     this->flat.T );

#endif
  }

  void Matrix_FLA_::solve_qr_2_r() {
#ifdef UHM_HIER_MATRIX_ENABLE
    // ----------------------------------------------------------
    // ** Hier-Matrix
    // ----------------------------------------------------------
    solve_qr_2_hier( this->fs, this->ss,
                     this->hier.ATL, this->hier.ATR,
                     this->hier.rt,  this->hier.rb);
#else
    // ----------------------------------------------------------
    // ** Flat-Matrix 
    // ----------------------------------------------------------
    solve_qr_2_flat( this->fs, this->ss,
                     this->flat.ATL, this->flat.ATR,
                     this->flat.rt,  this->flat.rb);
#endif
  }    


  static inline int solve_qr_1_flat( int fs, int ss,
                                     linal::Flat_ ATL, linal::Flat_ ABL, 
                                     linal::Flat_ t,   linal::Flat_ b,
                                     linal::Flat_ T ) {

    if (fs) {
      linal::Flat_ W;
      
      W.create    ( ATL.get_datatype(), T.get_m(), t.get_n() );
      
      FLA_Apply_Q_UT( FLA_LEFT, FLA_CONJ_TRANSPOSE,
                      FLA_FORWARD, FLA_COLUMNWISE,
                      ~ATL, ~T, ~W, ~t );
      
      W.free();
    }

    if (fs && ss) 
      FLA_Gemm( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE,
                FLA_MINUS_ONE, ~ABL, ~t,
                FLA_ONE, ~b );
    return true;
  }

  static inline int solve_qr_2_flat( int fs, int ss,
                                     linal::Flat_ ATL, linal::Flat_ ATR, 
                                     linal::Flat_ t,   linal::Flat_ b ) {

    if (ss && fs)
      FLA_Gemm( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE,
                FLA_MINUS_ONE, ~ATR, ~b,
                FLA_ONE, ~t );

    if (fs)
      FLA_Trsm( FLA_LEFT, FLA_UPPER_TRIANGULAR,
                FLA_NO_TRANSPOSE,
                FLA_NONUNIT_DIAG, FLA_ONE,
                ~ATL, ~t );

    return true;
  }

  static inline int solve_qr_1_hier( int fs, int ss,
                                     linal::Hier_ ATL, linal::Hier_ ABL, 
                                     linal::Hier_ t,   linal::Hier_ b,
                                     linal::Hier_ T ) {

    if (fs) {

#ifdef UHM_QR_INC_ENABLE
      {
        linal::Flat_ W;
        linal::Hier_ W_1;
        
        int bmn = get_hier_blocksize();
        
        W.create    ( ATL.get_datatype(), bmn, bmn*t.get_n() );
        W_1.create  ( W, bmn, bmn );
        
        linal::apply_q_inc( FLA_LEFT, FLA_CONJ_TRANSPOSE,
                                   FLA_FORWARD, FLA_COLUMNWISE,
                                   ATL, T, W_1, t );
        
        W_1.free();
        W.free();
      }
#else
      {
        linal::Flat_ W;
        linal::Hier_ W_1;
        
        int bmn = get_hier_blocksize();
        
        W.create    ( ATL.get_datatype(), T.flat().get_m(), t.flat().get_n() );
        W_1.create  ( W, bmn, bmn );
        
        linal::apply_q( FLA_LEFT, FLA_CONJ_TRANSPOSE,
			       FLA_FORWARD, FLA_COLUMNWISE,
			       ATL, T, W_1, t );
        
        W_1.free();
        W.free();
      }
#endif
    }

    if (fs && ss)
      linal::gemm( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE,
                          FLA_MINUS_ONE, ABL, t,
                          FLA_ONE, b);
    return true;
  }

  static inline int solve_qr_2_hier( int fs, int ss,
                                     linal::Hier_ ATL, linal::Hier_ ATR, 
                                     linal::Hier_ t,   linal::Hier_ b ) {

    if (ss && fs)
      linal::gemm( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE,
                          FLA_MINUS_ONE, ATR, b,
                          FLA_ONE, t );

    if (fs)
      linal::trsm( FLA_LEFT, FLA_UPPER_TRIANGULAR,
                          FLA_NO_TRANSPOSE,
                          FLA_NONUNIT_DIAG, FLA_ONE,
                          ATL, t );

    return true;
  }

}
